function [tval, prob, df] = tp_test(x,y)
%tp_test - Paired t test between two populations, assuming equal variance
% from Press et al., p. 618
% uses betai function
% Coded P. Manis, 12/18/2006


if(nargin == 0) % verify operation of the code with some data...
    x = [1:10];
    y = x;
    [tval prob df] = tp_test(x,y) % recursive...?
    y = x + 1;
    [tval prob df] = tp_test(x,y) % recursive...?
    return;

end;

tval = 0;
prob = 0;
df = 0;
if nargin ~= 2, % at this point...
    error('tp_test: Requires two input arguments');
end
% get rid of NaN's in input.
xn = find(~isnan(x));
yn = find(~isnan(y));
ok = intersect(xn, yn); % match where BOTH are not NaN
x=x(ok);
y=y(ok);

[m1 n1] = size(x);
[m2 n2] = size(y);
if (m1 ~= 1 & n1 ~= 1) | (m2 ~= 1 & n2 ~= 1)
    error('tp_test: Requires vector first and second inputs.');
end

n1 = length(x);
n2 = length(y);

if(n1 ~= n2)
    error('tp_test: Vectors must match in length for paired t-test');
end;

[mean_x, var_x] = mean_var(x);
[mean_y, var_y] = mean_var(y);
df = n1 - 1;

cov = 0;
for j = 1:n1
    cov = cov + (x(j)-mean_x)*(y(j)-mean_y);
end;
cov = cov/df;
svar = sqrt((var_x+var_y - 2*cov)/n1);
if(svar == 0.0)
    prob = 0.0;
    return;
end;
tval = (mean_x-mean_y)/svar;
prob = beta_i(0.5*df, 0.5, df/(df+tval*tval));


if(nargout == 0) % no output args, just print result
    fprintf(1, 't: %12.5f   p: %8.3f   df: %d\n', tval, prob, df);
end;
return
